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We studied formation of vortex with four-fold symmetry in a minimal model of self-propelled 
particles, confined inside a squared box, using computer simulations and also theoretical analysis. 

In addition to the vortex pattern, we observed five other phases in the system: homogeneous gaseous 
phase, band structures, moving clumps, moving clusters and vibrating rings. All six phases emerge 
from controlling strength of noise and contribution of repulsion and alignment interactions. We 
studied shape of the vortex and its symmetry in detail. The pattern shows exponential defect lines 
where incoming and outgoing flows of particles collide. We show that alignment and repulsion 
interactions between particles are necessary to form such patterns. Finally, we derived hydro- 
dynamical equations for our model and compared them with the results of both computer simulations 
and Quincke rotors. A good agreement between the three is observed. 
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I. INTRODUCTION 

All of us have seen fascinating movement patterns of 
flocks of birds m or schools of fishes aig. Simi¬ 
lar structures are widely seen in many places, ranging 
from our body size [6] down to nano meter scales mil], 
including either living individuals 0 HD] or non living 
ones mm- The common feature between all these 
diverse systems is activity among the individuals and 
therefore these systems are called active matter. Active 
matters, because of consumption and injection of energy, 
are always out of equilibrium and in recent years have 
attracted many attentions [TSl [14]. In one of the primary 
efforts, collective patterns from basic local interactions 
has been produced [15]. Later, phenomenological theo¬ 
ries as well as microscopic descriptions were established 
to find the characteristics and features of active mat¬ 
ter [T6H19] . In addition, more complex collective patterns 
were observed by introduction of new models [2Qll22] , and 
helped to improve knowledge of the phase transitions and 
behavior of active matter [231127] . 

Vortices are one of the most interesting patterns ob¬ 
served in active matter. Swirling of daphnia around a 
light shaft [28] , rotation of bacteria in droplets of bacillus 
subtilis morphotype which grow on an agar substrate [9] 
and vortices that are formed by the moving actin fila¬ 
ments on a surface coated with heavy meromyosin [8] 
are biological examples of vortex formation. There are 
also non-biological samples of vortex formation, e.g. ver¬ 
tically vibrated granular rods [29], anisotropic rods in 
a container [30] and micrometer sized insulator spheres 
well known as Quincke rotors m- In addition to stable 
vortex patterns, dissipation of vortices like a turbulent 
phase is observed in colonies of bacteria [32] [33] . Many 
studies have been done to understand vortex formation, 
its nature and characteristics. For example a model of 
self-propelled particles that repel each other in close dis- 
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tances and attract in far makes a giant vortex [34]. An¬ 
other example is the case of particles with intrinsic cur¬ 
vature in their motion like microtubules moving on a 
surface coated with myosines [35]. More sophisticated 
models consider chemotaxis and proliferation for bacte¬ 
ria to understand their swirling [9] . It is also possible to 
have a vortex array in the system [36] [37] , acquiring self- 
propelled particles with alignment and anti-alignment in¬ 
teractions with respect to distance [36] or time correlated 
noise [37] . 

In the experimental study of Quincke rotors [31] the 
particles and their interactions produce a complex vor¬ 
tex inside a squared box. The vortex shows a four-fold 
symmetry. This four-fold shape is along with an effect 
that we call hereafter “suppressed spreading”. In sup¬ 
pressed spreading, particles which are bounced back from 
a corner tend to spread over all available directions, but 
because of the flow of the other particles and collisions, 
the spreading of the outgoing flow is suppressed and is 
limited to a smaller angle. Suppressed spreading is visi¬ 
ble as a curved boundary, where the direction and density 
of particles change spontaneously, and we call it “defect 
line”. 

Inspired by the experimental vortex formations HHSH, 
here we study vortex pattern of self-propelled particles, 
confined in a geometry. Such a study has been done by 
simulation of active granular particles with inelastic in¬ 
teractions, confined in a squared box [38j|39]. However, 
these studies are limited in size and they could not de¬ 
scribe the complex structures and density jump lines of 
the Quincke rotors observed in the experiments [31] . 

First we introduce a minimal model and find the key 
elements required to have a vortex with the symmetry 
of confining geometry (section |^. The model is a simple 
generalization of continuous Vicsek model [9] with align¬ 
ment and repulsion interactions. Both interactions have 
physical interpretation and are derived theoretically in 
reference m- The important role of repulsion in vor¬ 
tex formation is revealed recently HOj. Next, we derive 
continuum hydrodynamic equations in the limits of high 
and low noise to compare the solutions with the parti- 
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cle model and experimental results (section [nl| ). Finally 
we present simulation results and discuss about patterns 
and their characteristic (IV). 


II. MODEL 


We consider two dimensional self-propelled particles 
with the same constant speed v. Angle of the velocity 
with X axis is 0 and the direction of motion for each 
particle is toward eg. The direction of each particle 
is changed by torque. This torque is originated from 
particle-particle and wall-particle interactions. For the 
dynamics of particles we consider 

ri=vee^, ( 1 ) 


Oi = Tf+T^+erii{t), (2) 


where fi is the position of the ith particle and eg. is a 
unit vector along swimming direction of the particle with 
angel Oi, {eg = cos{0i)ex + sin(^^)ey). In Eq. Q, rf and 
are the torques acting on the particle i from the other 
particles and the walls, respectively. We added a noise 
term er]i{t) which represents stochastic behavior of self- 
propelled particles and their environment, where r]i{t) is 
a Gaussian uncorrelated white noise with {r]i{t)) = 0 and 
{r]i{t)r]i{t')) = 6{t — t'), where e is the noise amplitude. 

The particle-particle interaction is a combination of 
alignment and repulsion. Alignment means particles ro¬ 
tate to make their moving directions parallel to each 
other and repulsion means that particles rotate to run 
away from each other (Fig.[^. The alignment and repul¬ 
sion interactions that are close to Quincke Rotors inter¬ 
actions m give. 


Tji = ^A{rji)x 
TT 

[(1 - a)sin {0j - 0i) + a{rji x • e^]. 


(3) 


The first term on the right hand side is aligning and the 
second term is repulsive torque, fji is the distance vector 
from particle j to i, Qp is the strength of particle-particle 
interaction and A{rji) is a function of interparticle dis¬ 
tance. 0 < (a < 1 controls the relative contribution of 
repulsion and alignment terms, i.e. a = 0 corresponds 
to the original continuous Vicsek model, and a = 1 cor¬ 
responds to a fully repulsive particle system. It is also 
interesting to see the result of negative a that is a combi¬ 
nation of alignment and attraction. The x sign between 
vectors shows vector product and . shows dot product. 
The model is two dimensional, but for a compact presen¬ 
tation we use dot product of z direction unit vector, Cz 
with the result of vector product as a scalar value. 

Figure shows a schematic presentation of interaction 
terms. As one can see, repulsion turns velocities of par¬ 
ticles in opposite of their inter distance direction. 



FIG. 1. (Color online) Repulsive torque between two particles 
in Eq. (|3 ) (a, b, c), and between a particle and a wall in Eq. § 
(d). Each particle is represented by a green disk. The straight 
red arrow shows direction of the particle. The arced orange 
arrow shows the direction of rotation, end of this arrow is 
particle’s final direction and its length has no relation to the 
magnitude of the torque. The thick vertical line is a wall. The 
dashed line in (a,b,c) is the inter distance line between two 
colliding particles and in (d) is the perpendicular line to the 
wall. In (a) particles are escaping and the torque intensifies 
this escape. In (b) particles are moving together and the 
torques cause them to move away. In (c) particles are moving 
toward one another and the torque inhibits them reaching 
each other. In (d) particle is moving toward the wall. 


Eor simplicity here we restrict ourselves to the case of 
constant a and consider a Heaviside step function for A 

A(rji) = e(R-rji), (4) 

where R indicates the range of interaction between the 
particles. 

Very similar to the repulsive interaction between par¬ 
ticles, if we label a wall by w, the applied torque by the 
wall w on a particle i is given by: 

'^wi ~ A{r \{finj X • ez \, (5) 

wi 

where is the unit normal vector of the wall, r^i is 
the distance of the particle from the wall, and is the 
strength of the particle-wall interaction. The factor l/r^i 
guaranties that the particles never pass the wall. 

One should note that in contrast to the alignment 
torque, the repulsion torque will not conserve the total 
angular momentum. In models with spontaneous vortex 
formation, a generating source of angular momentum is 
necessary. 

In the next section we will use this microscopic model 
to obtain hydrodynamic equations of the system. 


III. HYDRODYNAMIC EQUATIONS 

We can characterize the state of the system with den¬ 
sity of particles, p, and their polarization vector V = {eg) 
in space. The governing hydrodynamic equations of den¬ 
sity and polarization could be derived from microscopic 






equations. With the method presented in reference im, 
we derive Fokker-Plank equation for particle density up 
to second order of spatial derivatives, 
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df{0, f) 
dt 


= -(1 - a)gpde 


p27T 

f{0^ r) / dO' sm{0' — 0) 

Jo 




-de 

+ Dr 


agpR^ 

3 

56>2 


f{d,r)Wp{i^-{e0XeO 
vee ■ Vf{9,f}, 


(6) 


where f{0,f) is the density of particles at point Amov¬ 
ing in direction 6>, and = 6^/2 is rotational diffusion 
of particles. The terms on the right hand side represent 
alignment, repulsion, diffusion and advection of the par¬ 
ticles, respectively. 

In the following we solve the Fokker-Plank equation for 
two different regimes of high noise and low noise. 


A. High Noise Limit 

We use Fourier transform of density function /(6>, r) = 
where i is the imaginary number and fk 
is the kth Fourier component of /. Thanks to linear 
independency of Fourier basis we can split Eq. ^ to 
an infinite set of separate recurrence equations in Fourier 
space: 


dt 2 

{fk-iifi + 

-/,+i(/_i + :^v7-i)) 


agpR _ ^x(/fc_i - fk+i) 

^y{fk-l+fk+l)^ 

r\ fk — 1 f r\ fk — 1 f 

VOx -Z- VOy --- . 

2 ^ 2i 


( 7 ) 


Because the set of equations is unlimited we need to 
truncate it at some point. There is a damping term 
—Drk‘^fk with time scale Tk = —Drk‘^^ meaning that 
higher moments of fk vanish faster. Here we assume that 
moments of fk for k > 3 are zero and the second mo¬ 
ments converge to their equilibrium values fast enough 

that we can assume f ±2 = 0. This assumption is valid 
until the damping terms for higher moments are domi¬ 
nant. Comparing the coefficients of the right hand side 


of Eq. Q with ^ when /c = 3, we will get conditions. 
Dr ^ i^-oi)^9pR ^ ^ agpR ^ | to trun¬ 

cate Eq. 0 for k > 3. Given in the values of R = 1, 
a = O.b, Qp = 2 and v = 1 (the same as simulations) we 
find inequality Dr ^ 0.16 which satisfies all conditions. 

After the truncation, we can find f ±2 in terms of f±i 

and /o, then one can replace it in the equations of f±i. 
/o is density of particles and f±i is related to the po¬ 
larization of particles. Defining W = pV^ we can write 
f^i = Wx T 'iWy. Simplifying equations and solving the 
real and imaginary parts separately, it results in the con¬ 
tinuity equation, 

^+vV.W = 0, (8) 

and evolution of W, 


W = 


(1 - a)pgpE? 


- Dr 


(1 - afglR^W^ a^glR^ 

72Dr 

(1 - a)pgpR^ 


8Dr 


Vp 


w 


+ 


16Dr ' 16 

(1 - afglR^ 


V^W 


32Dr 


(w ■ v^tT) w 


_ (1 - a)gpvR^ r _ 5^t^2 

l6Dr ^ 2 

+ 5WV -W + SW- VW] 


+ 


a{l - a)glR^ 
12D^ 

appvR^ 


{W • \/p)W 


(2WV^p + 2Vp ■ VW 


48Dr 

- 3V ■ WVp + 3Vpx {V xW)]. 


( 9 ) 


On the right hand side of Eq. ([^, the first bracket con¬ 
tains driving terms which cause spontaneous polariza¬ 
tion, with the first three terms that usually appear in 
active matter hydrodynamic equations milSSl]. The 
fourth term introduces a reduction in polarity due to 
the net repulsive torque in density gradient. The sec¬ 
ond bracket is very important for spread of polarization. 
It is a diffusion-like term for W and comes from the align¬ 
ment of neighboring particle, closer than distance R. The 
third bracket comes from alignment interaction between 
particles and is well known in both phenomenological and 
analytical studies [161 SSHll] • Rest of the equation shows 
escape of particles from higher densities due to repulsion 
and advection. We are not worried about D~^ factor in 
the equation as we are in limit Dr ^ 0.16. 
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FIG. 2. (Color online) Results of continuum model [Eq. 
before occurrence of divergence with initial density po = 1, 
alignment interaction p = 2, repulsion factor a = 0.5 and 
noise Dr — 0.59 (e= 1.08) in a box with size of 120. |(a)| 
Density profile. |(b)| Velocity field, the lengths of arrows is 
proportional to the magnitude. Both images correspond to 
the same moment and their grid size is 32 by 32. 


Looking for steady homogeneous solutions of Eqs. 
and ^ with initial density po? can find a critical 
noise value Cc = {Dc = which 

the system behavior changes from polar to non-polar ho¬ 
mogeneous state. This change of behavior is also impor¬ 
tant for confined particles, because as simulations show, 
once the system goes to a polar state at this critical noise, 
particles start to rotate. With our simulation parameters 
one finds Dq = 0.5, therefore Eqs. and ^ are not 
able to explain the system for low noise values far from 
the transition point {Dr ^ 0.16). Solving the equations 
numerically, the answer for Dr > 0.6 is non-polar and 
homogeneous {W = 0) and with Dr = 0.59 numerical 
instabilities emerge and no vortex is observed. But inter¬ 
estingly, we can observe stripes forming and propagating 
in the box and reflecting from corners at the initial stage 
of computation before divergence occurs (Fig. |^. 

To obtain more stable hydrodynamic equations we de¬ 
rive them in the low noise limit. In the low noise regime, 
density and polarization change slowly in space, that en¬ 
hance the stability of equations. 


(6(6>)), we find. 


X {b'{6)sm{6' - 9))0>fi 

X + ^p(f)VV(^ 

+ 1(1 - a)gR'^p{rf‘ 

X {h'{e)y\sm{e' -6))0,)0 

+ 1(1 - a)5i?Vr)W(rl 
■{h'{e)V{sin{e' -9))e')0 
^9pR 


( 10 ) 


+ 


p(f)Vp(f) . {{b\0)eo) X e,) 


-v\/-{p{r){b{ 0 )eo)) 

+ ^(6(^))V-(p(f)(e,)) 
^Drp{r){b"{0))+O{\/^), 


where indicates average over 0 with probability dis¬ 
tribution P{0^f). In the low noise limit, P{0^r) is 
sharply peaked around the mean and we can approxi¬ 
mate sm{ 0 ' — 0 ) with 0 ' — 0 to find homogeneous solu¬ 
tions of Eq. (10) for b{0) = 0 and b{0) = This gives 
us dynamical equation for variance of 


-^ = Dr - {1 - a)gpPocrj, ( 11 ) 


B. Low Noise Limit 


In this limit, we decompose f{0,f) to the number den¬ 
sity p(r) and orientational probability distribution of par¬ 
ticles, P{0^f) [i.e. f{0,f) = P{0,r)p{f)]. We multiply 
both sides of Eq. by an arbitrary function b{ 0 )^ and 
integrate over 0 to find the governing dynamics of {b{ 0 )). 
Then for 6(6>) = 1 it gives us continuity equation in terms 
of {ee) = V and p [p = —vW • {p{r){e 0 ))]. Substituting p 
from continuity equation into the governing equation of 


where Dr is acting here as a source of dispersion for 6>, 
while alignment interaction (1 — a)gpp reduces the dis¬ 
persion. From this equation we find that decays with 
time scale r^r = 1/(1 — OL)gpp to its equilibrium value 
(7^ = Dr/{l—a)gpp. Our assumption for sharpness of dis¬ 
tribution function is true when <C 1 and this gives us a 
limit for the noise. Dr <C {l — a)gpp. By expanding up 
to second order around the mean value 0 ^ and using the 
facts that = Px '^PPy and = ^{Vx PP-PPy), 

we can find a relation between and polarization P as, 
P = {1 — (7^ /2). The same assumption - small deviation 
- in Eq. (10) with b{0) = helps to find the dynamics 
of the polarization, P{rR), 
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dt 


= V 


2(1 - a)gR^ 


y.{p{^ + ^V^p{m^-V)-D, 

+ 1(1 - a)gR^p{r) \{2V - l)V‘^V 


- (4P - 3)PP ■ 


+ -(1 - a)gR^ [(2P - l)Vp(f) ■ VP 


4P-3 

2P2 

agpR^ 


PVp{f) ■ Vp2 

(1 - 2P)Vp 


( 12 ) 


(4P - 3)P • V/op] 


+ V 

1 


-V(p(P-l)) 


+ ^P ■ V (p(p - 1)(P - 3)p) 

+ (P - 1)(P - 3)PV ■ P - P ■ VP 
+ 0(V3). 


The first bracket on the right hand side of Eq. (12) pre¬ 


vents V to become zero, the second bracket spreads po¬ 
larization because of interaction with neighboring parti¬ 
cles, the third bracket shows an alignment competition 
between high and low density regions, the fourth bracket 
represents repulsion of particles (moving against Vp) and 
the last one shows advection. 

It is clear that solving Eq. together with continu¬ 
ity equation in squared geometry is complicated, there¬ 
fore we restrict our calculation to circular box where the 
scalar variables depend only on radial component r. We 
also neglect second order derivations that are related to 
shorter changes of fields. Even with these assumptions, 
presence of a non-homogeneous steady state is not ob¬ 
vious and we need to write everything up to the first 
order of Dr/{1 — ot)gpp- With all these simplifications 
one finds that continuity equation implies that polariza¬ 
tion is toward polar coordinate unit vector 0, and the 
density satisfies the following relation: 


/X Dr , . . 3v f r 

+ n -^ = — In — 

{l-a)gp ag \ro 


(13) 


where ro is a length scale related to the initial density 
of particles. We can see that at distances shorter than 
ro the density gets very small values. Similar results 
with logarithmic dependence of density is obtained in 
Quincke rotors m- If the system has no noise {Dr = 0) 
the right hand side of Eq. (13) is negative for r < ro that 


corresponds to a zero density in center. In fact no particle 
reaches the center in noiseless system, but when turning 
on the noise some particles could reach the center with a 


small chance. These results are in good agreement with 
simulation results (see Eigs. |^and[^in section IV). 

Integrating equation © for Dr = 0 over the box and 
equating it with the total number of particles, one obtains 
an expression for ro, 


Po = 


3v 

2(xgp 


p2 

v^box 


- 1 - 2 In 


_ro_ 

-^boi 


(14) 


where Rbox is the radius of the circular box. Eq. (14) 
shows ro decreases either by increasing po or agp, both 
causing stronger steric repulsions. Although our dynam¬ 
ics is not Newtonian, we see dependence of ro/Rbox on 
V. Having a mixture of particles with different velocities, 
this effect can separate fast and slow particles. We will 
discuss this effect and phase separation of slow and fast 
particles in a future work. It is interesting to mention 
that a mixture of fast and slow particles could be pre¬ 
pared experimentally by using different sizes of Quincke 
rotors due to their velocity relation. 


j{Eo,a)=aC^{Eo/EQP 


- 1 , 


(15) 


where v is the velocity of particles, a is their radius, C 
is a combination of elements of modified mobility matrix 
and Maxwell-Wagner time of the environment, Eq is the 
applied electric field, Eq is a threshold for electric field. 
In this relation Eq and C depend on both environment 
and a 1^. 


C. Numerical Method in Computation of Noiseless 
System 

To find steady state solution of Eq. ( [T^ , when Dr = 0 
{V = 1), we need to rewrite Eq. (12) i n term s of W to 
avoid numerical instabilities (subsection [III A ), 


W=-{l-a)gpR^ 

X {piV'^W 

-1—K^)(e 2 xp){e^xp)- Vp{r) 

-v{p^ -W + W 0{V^). 

To compare theoretical result with simulations and ex¬ 
periment, we also numerically integrate Eqs. (16) and 
(§ with a square box geometry. Integration of the 
equations are done using pseudo spectral method, semi- 
implicit time stepping and anti-aliasing (2/3 rule) tech¬ 
niques [45j |46]. Using these techniques alone, is not 
enough to have stability and we need to add a p 
{K = 0.2) term to the continuity equation to avoid nega¬ 
tive density values. This extra term does not drastically 
change the underlying physics of the model and this tech¬ 
nique was used before in the literature m- Additionally 
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FIG. 3. (Color online) Density of particles p{r) as a function 
of radial distance r, for a noiseless system in a circular box 
with diameter L = 120 and initial density po = 8. Labels show 
different values of a. One can see that the density increases 
logarithmically with r [Eq. |T^]. 

at the initial stage of integration we must set v = 0.005 
and slowly increase the speed up to the highest stable 
value which is v = 0.3. We set a = 0.5 and Qp = 2 
to compare our result with the particle model. A circu¬ 
lar rotating pattern with homogeneous density for initial 
condition is necessary for stability. With the use of image 
method we apply slip boundary condition. That means 
we replicate our system on the sides with reflecting W, 
to justify slip boundary condition. 

Now we have the right theoretical tools to compare our 
simulation results. In the next section we will provide 
the method and the results of simulations based on the 
microscopic model to compare with the theory. 


(a) 



(b) 



IV. SIMULATION 

To simulate our model, we used integration technique 
(Ito method [47]) with time steps dt = 5 x 10“^. We set 
R = 1, Pp = 2, =40. The speed of particles is set 

equal to one {v = 1), unless otherwise stated. The con¬ 
trol parameters in simulations are the strength of noise 
e, repulsion strength a, and initial density po- Except 
for part |IV A] where we use a circular box, we always put 
N particles inside a squared box with four surrounding 
walls of size L. Simulation box must be large enough to 
capture all aspects of a self-propelled system. This limit 
is originated from band structures j48U5Q] that have typ¬ 
ical length scale of 'T’/e^. Then for given parameters, one 
finds the condition L ^ v/cc ^ 4. We initially posi¬ 
tion particles homogeneously in a triangular lattice with 
uniformly random direction of motions. We run the sys¬ 
tems long enough to guarantee that the effects of initial 
condition vanish, and then sample the systems. 


FIG. 4. (Golor online) Gomparison of simulations (points with 
error bars) with theoretical Eqs. (13) and (14) (dashed lines) 
of a circular box with radius Rbox = 60 and different initial 


densities. |(a)| The plot shows slope in Eq. (13) as a function 
of a. |(b)|Th e plot shows ro as a function of a with different 
po [Eg. (l4[)]. 


A. Circular Box 


First we use a circular boundary to compare with the 
solutions of our theoretical Eqs. (13) and ( [l4| ). Figure]^ 
shows particle density as a function of radial coordinate, 
r, for noiseless microscopic model in a circular box with 
radius Rbox = 60^ and initial density po = 8. The fig¬ 
ure shows the results for different values of a. Density 
is increasing logarithmically with r, and by increasing 
Of, reduction in the slopes and the size of empty region 
in center is visible. Equation (13) predicts exactly the 
same behavior. For a more precise comparison we also 
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FIG. 5. (Color online) Snapshots of the simulations with var¬ 
ious noises and fixed repulsion factor a = 0.5, initial density 
po = 1 and box size of L = 120. Open circles indicate par¬ 
ticles, color shows the direction of motion corresponding to 
the color wheel located at top left corner of each box, and 
particles tails are along their trajectories. Squared windows 
represent zooming in a part of the box, the length of selected 
area and the magnification factor is written close to the cor- 
responding wi ndow . By decreasing the noise, one can see 
homogeneous, [(1^ inhomogeneous with curved stripes and 
vortex. In the vortex phase, suppressed spreading and four 
defect lines are observable. 


plot the slopes and vq values versus a for given density 
values p = 1, 2,4, 8 in Fig. Dashed lines in Fig. show 
Eqs. (13) and (14) for each given density. We observe a 


good agreement between theory and simulation specially 
for higher densities. This is because in high density our 
assumption of having sharp distribution for 0 is more 
accurate. 


B. Phases 

During simulations by changing control parameters e 
and a, we observed six different regimes (phases): homo¬ 
geneous gaseous, band structures, moving clumps, mov¬ 
ing clusters, vibrating rings and vortex. At a given 
value of a = 0.5, and by decreasing noise, homogeneous 
gaseous phase, band structures and vortex pattern are 
observed (see Fig. [^. Starting fro m very high noise, the 


system is in a gaseous phase [Fig. ^[a) . By decreasing 


FIG. 6. (Golor online) Snapshots of the simulated system with 
po = 1, e = 0.1 and different values of a. Dots indicate parti¬ 
cles, color shows the direction of motion corresponding to the 
color wheel located at top left corner of each box. Squared 
windows represent zooming in a part of the box, the length of 
selected area and the mag nifica tion factor is written close to 
the corresponding window, (a) The system with a very small 
value of a where particles rotate cohesi vely around the box 
and a very small defect line is observable. |(b) [ Switching off the 
repulsion (a = 0) a clump forms with large number of parti¬ 
cles in a small area. This clump is moving and gets reflected 
from the walls while at some times it is divided to smaller 
clumps due to the noise a nd a t other times small clumps join 
to make a bigger clump. |(c)| Setting a = —0.01 cluster s of 
particles continuously divide and join during simulation. |(d)| 
By further decrease in a vibrating rings emerge. 


the noise strength, the density gets inhomogeneous and 
particles form traveling curved stripes which are due to 


the reflections from the walls and the corners [Fig. m 
By further reducing the no ise, p articles start to rotate 
in the box as shown in Fig. The direction of rota¬ 

tion is random and depends implicitly on initial positions 
and the string of random number samples. In rotation 
[Fig-IE], one can easily recognize suppressed spreading 
and the presence of the defect lines. The same pattern 
has been experimentally observed in the suspension of 
Quincke rotors m- 

Moving clumps, moving clusters and vibrating rings 
could be observed by changing the ratio of repulsion to 
alignment (Fig. [^. For this purpose we set e = 0.1 and 
change a. If we have no alignment {a = 1), we obtain 
a gaseous homogeneous state again, but particles in this 
homogeneous state have a more robust ballistic motion 
in comparison to the homogeneous state observed in high 
noise. A slight de crease in a could lead to formation of 
a vortex, e.g. Fig. a) [ shows the rotation for a = 0.01. 
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FIG. 7. (Color online) Velocity autocorrelation over time for 
systems with a — 0.5, and various e. (a) The noise 


is lowe r tha n critical noise and the system is in the vortex 
phase. |(b)| A big polar structure rotating (e = 0.6), locally 
propagating band structures (e = 0.7, 0.8,1.0), and homoge¬ 
neous gaseous phase (e = 1.2). 


Since the repulsion between particles is not strong, they 
form dense bunches near the walls with a large empty 
space in the center in such a way that almost all the 
particles walk on the walls. Increasing the repulsion, 
more space is covered by the particles. By turning off 
the repulsion completely, particles form very high den¬ 
sity clumps that bounce off the walls [Fig.[f 


Clumps 

formation is not observed in periodic boundary condition, 
because the particles do not meet each other as frequently 
as confined particles. The shape of clumps in a squared 
and circular box is square and circle respectively. For 
negative a’s, particles rotate toward each other, and if 
this attraction is small we see multip le clusters traveling 
and bouncing off the walls [Fig. m- These clusters are 
extremely packed and unstable. They may divide into 
smaller groups or join together to make a bigger mass 
of particles. The division takes place when particles are 
more distant than their interaction range so that they 
cannot return to each other. Finally a strong attraction 
produces vibrating rings of particles when particles can 
not escape from the ring [Fig. 11. 


To find detailed properties of these phases we look at 
velocity autocorrelations over time (C^(t) = {vi{t)'Vi{t^ 
and space {Cv{r)). We see that autocorrelations 
in time oscillate in both vortex phase and band struc¬ 
tures (Fig. [^. The oscillation in vortex phase [Fig. m 
is damping over time gradually. This damp is due to 
the different frequency of rotation of particles. For in¬ 
stance in a vortex in a circular geometry, Cy(r) is com¬ 
puted by integrating velocity autocorrelation of individ¬ 
ual particles: Cv(r) = 1/N 27rrp{r) cos{vr/r)dr^ 

where rmin and Vmax are minimum and maximum radii 
of rotation respectively and v is the speed of particles. 
One can see that in large r, the cosine in the integral 
changes fast and if p{r) is a smooth and slow changing 
function, positive and negative parts of the integral in 
one oscillation cancel each other. Hence, the integration 
gives us a negligible result when r is large. Instead, in 
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FIG. 8. (Golor online) Velocity autocorrelation over dist ance 
for systems with po = 1, a = 0.5, and various e. |(a)| The 
noise is lower than critical noise and the system is in the 
vortex phase. |(b)| A big polar structure rotating (e = 0.6), 
locally propagating band structures (e = 0.7,0.8,1.0), and 
homogeneous gaseous phase (e = 1.2). 



0.9 1 


' p=0.2 ' 

+ 


0.8 ' 

■ □ 

p=0.4 

X 

- 

0.7 

X 

^ Sx 

p=0.6 

p=0.8 

□ 

- 

0.6 

□ n 

+ 

p=l 

■ 

- 

0.5 

□X + 

□ 



- 

0.4 

□ X 



- 

0.3 

■ ■ 



_ 


X + 




0.2 

■ 



- 

■ 




0.1 

X X + 



- 


Dl m — a ^ ^ ^ ; 


0 0.2 0.4 0.6 0.8 1 1.2 1.4 

e/Ec 

FIG. 9. (Golor online) Scaled angular momentum per particle 
(2M/L) versus scaled noise e/cc for different values of density 
pQ. Simulation box size is L = 120. By decreasing noise, 
system shows a phase transition. Moreover, data for different 
densities collapse on the same curve. 


the traveling band st ructur e phase, C^(r) shows more ro¬ 
bust oscillations [Fig. 7[[b) . Band structures travel in the 
medium with a constant speed which depends on noise 
strength. These waves bounce from the walls and pe¬ 
riodically move from one corner to another. In homo¬ 
geneous gaseous state, Cv{r) spontaneously reaches zero 
[Fig.[lb)1e=1.2]. 

In vortex phase C^(r) (Fig. shows the behav ior and 
characteristic length scale of rotation [Fig. m- This 
length scale is the same as the box dimension. In higher 
noise no negative correlation in long range dis tance s is 
seen except for rotating pattern in e = 0.6 [Fig. m 


We claim that the vortex formation is a consequence 
of the confinement of polar state. This claim is sup¬ 
ported by the coincidence of the vortex transition point 
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FIG. 10. (Color online) Phase diagram of the system with 
Pq — 1 and L = 60. This diagram is constructed from the 
results of 625 points, each averaged over 15 realizations. All 
boundaries are constructed by looking at angular momentum 
and local cohesion except the boundary between rotation and 
disordered phase which is from theoretical prediction. Rota¬ 
tion is recognized by looking at angular momentum per par¬ 
ticle greater than 0.5. Vibrating rings where recognized by 
local cohesion smaller than 0.1 when a < 0. The remaining 
points are either in the moving clusters or the clumps phase. 
All points with a = 0 correspond to formation of moving 
clumps and the rest of points are in moving clusters regime. 


of confined particles with the non-polar to polar transi¬ 
tion point of unconfined particles. To characterize vortex 
transition we use the average angular momentum of the 
particles, 



N 

vRi X 60, 


(17) 


Since M scales with L, we plot 2MIL as a function of 
scaled noise ejcc in Fig. We observe a transition in 
angular momentum at e = Cc for all densities. 

The phase diagram of the system with six different 
phases is shown in Fig. To define vortex phase we 
looked at angular momentum per particle ((M) > 0.5). 
To recognize vibrating rings and disordered phases, we 
computed local cohesion. 


$ = 




^ cos(6>i - 9j) 


(18) 


Here is the number of particles in a certain window 
and the sum is over particles within that window. 



FIG. 11. (Golor online) Time average of p, $ and v in space 
for po — 1 and e = 0,0.5 in a box of size L — 120. |(a 


and 1(b) I density with grid dimension 128 by 128. |(c)| and [(d 
cohesion with grid dimension 64 by 64. (e) and |(f)| velocity 
with grid dimension 32 by 32. Golor codes corres pond to the 
magnitude of each held a nd the le ngth of a rrows in |(e)| and |(f)| 
shows the average speed, (a H and 1(e) I belong to the s ame 
nois eles s simulation of a clockwise vortex (e = 0). (b)[ (d)[ 
and 1(f) I belong to the simulation of a counter clockwise vortex 
with noise (e = 0.5). Defect lines corresponding to suppressed 
spreading behavior are observed in reagion with low cohesion 
and high density. 


C. Suppressed Spreading and Defect Lines 

Suppressed spreading and presence of the defect lines 
are two of the interesting features of the system. For a 
clearer observation of defect lines we could measure the 
cohesion between particles. Figure shows time aver¬ 
aged p(r), ^(f^, and v{r). Near the corners, density gets 
its highest value. Velocity held shows an outgoing how 
of particles at each corner which is suppressed by colli- 
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(a)particle simulations (b)continuum model 


(a)p, <l> along yo = 53 


(h)v along yo = 53 



FI G. 1 2. (Color online) C omp arison of average density field 


(c)a; on line yo = 0 


(d)v on line yo = 0 


for (a) particle model and |(b)| continuum model. Parameters 
are set to po = 1, S' = 2, a = 0.5, e = 0, and L = 120. The 
color represents density. The grid size is 128 by 128. Few 
points of continuum model are in range p G (—0.1,0) but 
to have a better comparison we plot both maps in the same 
interval p G (0, 2.5). 


(a)particle simulations 


(b)continuum model 


0.4 
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FI C. 13 . (Color online) C ompa rison of average velocity field 
for |(a)| particle model and |(b)| continuum model. The color 
represents velocity magnitude. The length of arrows is pro¬ 
portional to the speed. The grid size is 32 by 32. To enhance 
the velocity map of the continuum model, we removed points 
with density lower than 0.03 that have large computational 
errors. 


sion of incoming flow to the same corner. This collision 
makes a defect line that corresponds to a lower cohesion 
but high density. The name defect line is used, due to the 
spontaneous change of velocity in the boundary between 
incoming and outgoing flows. 

Defect lines are also present in the solution of the con¬ 
tinuum model introduced in subsection III C| Density and 
velocity fields in the numerical solution and particle sim¬ 
ulations are plotted in Figs. and respectively, for 
a system with po = 1? = 2, a = 0.5 and v = 0.3 

which show excellent qualitative but not precise quanti¬ 
tative agreement between the simulation and the contin¬ 
uum model. 

Figure shows the change in average density, cohe¬ 
sion, velocity and rotational frequency of particles across 
a horizontal line y = yo for yo = 53 and yo = 0. One 
can see that right after defect (lower cohesion), density 


3 

2.5 

2 

3 1.5 
1 

0.5 
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-0.4 -0.2 0 

x/L 


0.2 0.4 



FIG. 1 4. ( Color online) Cross sectional value of ti me a vera ged 
fields, (a) p and 4> across horiz onta l line y = 53. (b)| and |(d) 
Vx and Vy across horizontal lines |(b)| y = 53, and |(d)| y — 0. |(c) 
Angular velocity u across horizontal line y — d. Simulation 
parameters are po = 1.0, e = 0.05 and L = 120. Sub-caption 
of each plot shows the value of yo and the corresponding value. 
One can see that density exponentially decreases and the co¬ 
hesion has a fall-off corresponding to crossing the defect line 
Velocity component perpendicular to the wall is con- 


(a) 


stant along th e wa ll [(b)[ an d it l inearly increases with distance 
from the wall |(d)| Finally (c) shows that rotation frequency 
of particles depends on their position. 


is exponentially decreasing [Fig. l^[a)| . By looking at the 
velocities we find that velocities are tangent to the wall 
and by getting away from the wall we see that perpen¬ 
dicular component of veloc ity to the wall is independent 
of horizontal position [Fig. fT3Fb)Fand is proportional to 
the distance from the wall [Fig. l^[d)] . This proportion¬ 
ality corresponds to the exponential shape of the defect 
line. To show this, we first defi ne u{x ) the distance be¬ 
tween the top defect line in Fig. from the top wall 

as a function of x. Linear relation of perpendicular com¬ 
ponent of velocity to the wall gives us Vy = au where a 
is a positive number. This linear relation lets us find a 
differential equation with exponential solution: 


du{x) 

dx 


= —= 


1-vl 


—Vy = au{x). (19) 


Here we used t he app roximation that \vy\ 1 which is 
clear from Fig. m\ _ 

Defect lines in the experiment of reference m also 
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x/L 


FIG. 15. (Color online) The defect line distance from the 
wall 1 / as a function of horizontal position x from the corner 
in three cases of Quincke rotor experiment, simulation and 
continuum theory. Simulation and continuum parameters are 
set to po = 1, = 0.5, V — 0.3,e = 0, and L — 120. Experi¬ 

mental data are averaged movie frames of reference m with 
enhancement of final image to avoid computational error and 
whiteness considered proportional to density, u is dehned as 


the points where - ^ >0.1. One can see that at 

\P ^y\(x,L/2-u{x)) 

the interval —0.3 < x/L < 0.2, the defect line shape is expo¬ 
nential which is in agreement with the exponential decay of 
density in Fig. 


We also see that experiment, simulation 


and continuum model are all in good agreement. 


seem to have exponential form. We extracted experi¬ 
mental images from a movie in supplementary material 
of Quincke rotor experiment [31] and assumed the aver¬ 
aged gray scale over the images is proportional to the 
density. The form of defect lines in simulation, contin¬ 
uum model, and experiment are sketched in Fig. that 
shows a very good agreement between these three data 
sets. _ 

Figure l^[c)| shows that the angular velocity of par¬ 
ticles across y = 0 is not constant. We also computed 
the average angular velocity {uj = AO/At) of any parti¬ 
cle around the center of the box and observed that the 
angular velocity is almost proportional to the inverse of 
distance from the center {uj ^ 1/^)- That means parti¬ 
cles velocity component perpendicular to their position 
vector is rather constant {uj = v±/r). The result shows 
that in the low noise regime v_\_ ^ {v) ^ v = 1; However, 
a reduction in v± as well as the time averaged velocity of 
the particles is observed when the noise is high. 


V. DISCUSSION 

In this study we presented a minimal model to mimic 
the behavior of Quincke rotors in a squared box. Our 


model has alignment and repulsion between particles and 
it shows six different phases that are homogeneous disor¬ 
dered, moving clumps, moving clusters, vibrating rings, 
traveling stripes and vortex. Our focus in this paper 
was to study the vortex formation. We derived hydrody¬ 
namic equations in high noise and low noise limits and 
saw that each of them could describe some aspect of the 
system. With high noise equations we are able to predict 
the transition point while we are not able to find a vortex 
formation; However, low noise equations give us vortex 
solutions. The shape of vortex resulting from continuum 
model is very similar to our simulations and the Quincke 
rotors experiment. In the theoretical model, simulations 
and experiment, the center of box is rather empty and we 
observe a four-fold symmetry in the shape of the vortex. 
This four-fold symmetry is because the spread of particle 
flows going out of a corner is suppressed by collision of 
another flow and a defect line emerges for each corner. 
The shape of this defect line is exponential in experiment, 
theory and simulations and we brought some evidences to 
prove it. In addition to the shape of defect line, our the¬ 
oretical calculation well predicts the size of empty region 
in the center of the box and shows that the size depends 
on repulsion, total density of particles and their veloc¬ 
ity. The velocity dependence is important and one could 
use self-propelled particles inside a box, as separator of 
fast and slow moving particles. Because Quincke rotors 
velocities (Eq. depend on their radius and other envi¬ 
ronmental factors m, they are particular candidate for 
future experiment about this separation technique. Our 
next aim of research is to study the behavior of mixture 
of fast and slow particles, using simulations. Other stud¬ 
ies of similar systems could be done like a more precise 
study of transition with finite size scaling or a study of 
other geometries of the boundary, e.g. n-sided polygon. 
One could also find the limit of n in an n-sided polygon in 
which the vortex does not have the symmetry of the poly¬ 
gon anymore. Finally, one can construct hydrodynamic 
equations using Gaussian approximation and check the 
accuracy of the results m- 
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